
%***********************************************************************
% wtd_moms.m
%***********************************************************************
% C.-Y. Cynthia Lin, 12 February 2010
%-----------------------------------------------------------------------
% The code provided is from the following paper:
%
% Lin, C.-Y. Cynthia.  (2013).  Strategic decision-making with information 
% and extraction externalities:  A structural model of the multi-stage 
% investment timing game in offshore petroleum production.  Review of 
% Economics and Statistics, 95 (5): 1601-1621.
%
% If the code is helpful to you in writing your own code, I would appreciate
% it if you would cite the above paper in your work and acknowledge me 
% for help with the code.  Thank you!
% ------------------------------------------------------------------------
%
% This function calculates the weighted moment condition for the
% oil production model given the specified choice of estimators for Vce.
% # moments = # params
%
% -----------------
% MOMENT CONDITIONS
% -----------------
% The moment function for each market-time observation j is: 
%    f(data_j,params) = ...
%      [ ge(tuple_j, params) - ge_emp_j(tuple_j) 
%	 gd(tuple_j, params) - gd_emp_j(tuple_j) 
%	omega_j * ( ge(tuple_j, params) - ge_emp_j(tuple_j) )
%	omega_j * ( gd(tuple_j, params) - gd_emp_j(tuple_j) )  ]
%
%
% The population moment conditions are given by:
%     E[f(data_j,params)]=0 
%
% The sample moment conditions are given by:
%    (1/n sum_j{f(data_j,params)})=0 
% and amount to:
%    sample_mom_conds = 
%    (1/n sum_j{f(data_j,params)})=
%      1/n *
%       [ sum_s{tot_not_e_tuple(s) * (ge(s) - ge_emp(s))} 
%         sum_s{tot_e_not_d_tuple(s) * (gd(s) - gd_emp(s))} 
%         etc     ]
% where s indexes the tuples reached in the data
%
% The weighted moment conditions is (1 X 1)
%    (sample_mom_conds)' weight_matrix^(-1) (sample_mom_conds]) 
% where
%    weight_matrix =  1/n sum_j{f(data_j,params)f(data_j,params)'}
%
%
%
% =========================================================================
% ----------------
% Input Arguments:
% ----------------
% 	params
%		vector of parameters to be estimated (k X 1)
%		   params =
%			sigma_mu
%		 	sigma_eps
%			ce_factor
%			gamma_tote
%			gamma_totd
%			gamma_bid_bin
%			gamma_cd_bin
%			gamma_oilp_bin
%			gamma_constant
%	lease_term
%		lease term
%	Te
% 		Last period of consideration for the exploration stage
% 		(this is the last period in which decisions can be made)
%	Td
% 		Last period of consideration for the development stage
% 		(all years from lease_term onwards are lumped together
% 		b/c tot_e_t, which depends on t, will not change from then
% 		onwards)
%	max_mkt_size
%		maximum market size (scalar)
%	beta
%		discount factor
%	tuple_mat
%	   	matrix of all the tuples of state variables that appear
%	   	in the data 
%	   (N_TUPLES X N_STATE_VARS)
%	n_per_tuple_vec
%	   	number of times ecah tuple of state variables was
%		reached in the data
%	   (N_TUPLES X 1)
%       tot_not_e_vec
%	   = total # who have not explored
%	     before year t when the state tuple at time t
%            is "tuple"	
%	   (N_TUPLES X Te_m)
%       tot_e_not_d_vec
%          = total # who are explored at t
%            (ie, by (t+1)) 
%   	     but have not developed
%	     before year t when the state tuple at time t
%	     is "tuple"	
%	   (N_TUPLES X Td_m)
%	n_d_vec
%	  = # who develop at time t 
%	    (ie, between t and (t+1)) when the state tuple 
%            at time t is "tuple" 
%	   (N_TUPLES X Td_m)
% 	Mn_emp
% 	 	empirical transition matrix from the point of 
% 		view of the owner of an unexplored tract 
%		that decides not to explore at 
% 		at time t.  The (i,j,t)th element gives the 
% 		empirical probability that the state
% 		tuple next period is tuple j, 
%		given that the tuple this period is
% 		tuple i and given that this tract owner 
%		decides not to explore this 
% 		period.
%	   (N_TUPLES X N_TUPLES X Te_m-1)	
%       Me_emp
% 		empirical transition matrix from the point of 
% 		view of the owner of an explored tract 
%		that decides not to develop at 
% 		at time t.  The (i,j,t)th element gives the 
% 		empirical probability that the state
% 		tuple next period is tuple j, 
%		given that the tuple this period is
% 		tuple i and given that this tract owner 
%		decides not to develop this 
% 		period.
%	   (N_TUPLES X N_TUPLES X Td_m)	
%       ge_emp
%		empirical probability of exploring 
% 		at time t given that the tract
% 		has not yet been explored and 
%		conditional on the publically
% 		available information at time t 
%		(but excluding the idiosyncratic
% 		private shock eps_i_t).
%	   (N_TUPLES X Te_m)	
%       gd_emp
%		empirical probability of developing
% 		at time t given that the tract
% 		has been explored but not yet developed 
%		and conditional on the publically
% 		available information at time t 
%		(but excluding the idiosyncratic
% 		private shock eps_i_t).
%	   (N_TUPLES X Td_m)	
%	avg_profitsifd_emp
%		average gross profits conditional on developing
%   		for each tuple reached in data
%	   (N_TUPLES X 1)
%	weight_matrix
%		matrix of weights for the moment conditions
%	   (2*N_TUPLES X 2*N_TUPLES)
%       tol_fp              
%		termination tolerance when determining the
%		fixed point for value functions 
%	ini_diff_fp
% 		Initial difference between value functions before 
% 		solve for fixed point
%	Vce_estimator
%		one of:  'I' or 'G'
%
% -----------------
% Output Arguments:
% -----------------
%       wtd_moments             weighted moment conditions (1 X 1)
%
%
%
%**************************************************************************
function wtd_moments = ...
    wtd_moms(params, ...
    		lease_term, Te, Td, ...
    		max_mkt_size, ...
		beta, tuple_mat, n_per_tuple_vec, ...
                tot_not_e_vec, tot_e_not_d_vec, ...
		n_d_vec, ...
          	Mn_emp, Me_emp, ge_emp, gd_emp, ...
		avg_profitsifd_emp, ...
   		weight_matrix, tol_fp, ini_diff_fp, ...
		Vce_estimator);

% Global variables
    % Estimated continuation values
    global Vcn
    global Vce

    % Estimated investment probabilities
    global ge
    global gd

    % Estimated deterministic component of profits
    global pi0_e
    global pi0_d

    % Estimated ex ante values
    global EVn
    global EVe

    % Weighted moments
    global wtd_moments



% Matlab index of last period of consideration for each stage
% (this is also the length of the time dimension in each stage)
Te_m = Te + 1; 
Td_m = Td + 1; 

% Determine n, the number of observations.
n = sum(n_per_tuple_vec); 

% Determine k, the number of parameters.
k = length(params);

% Determine the number of different tuples reached in the data
n_tuples = size(tuple_mat, 1);

% Determine the maximum number of tracts in any tuple
max_n_per_tuple = max(n_per_tuple_vec);


% Separate params into its individual elements    
sigma_mu	= params(1,1);
sigma_eps	= params(2,1);
ce_factor	= params(3,1);
gamma_tote	= params(4,1);
gamma_totd	= params(5,1);
gamma_bid_bin   = params(6,1);
gamma_cd_bin  	= params(7,1);
gamma_oilp_bin  = params(8,1);
gamma_constant  = params(9,1);


% Separate tuple_mat into vectors for each state variable
tot_e_t_vec		= tuple_mat(:, 1);
tot_d_t_vec		= tuple_mat(:, 2);
bid_bin_t_vec		= tuple_mat(:, 3);
ce_bin_t_vec		= tuple_mat(:, 4);
oilp_bin_t_vec		= tuple_mat(:, 5);

% Set cd_bin = ce_bin
cd_bin_t_vec = ce_bin_t_vec;


% If sigma_mu or sigma_eps are not strictly positive, then
% set the weighted moment conditions to be equal to infinity
% and return to the calling program
if (sigma_mu <= 0) | (sigma_eps <= 0)
    wtd_moments = Inf;
    return
    disp('Error: did not return');
    break
end

% If any of the params are NaN, then 
% set the weighted moment conditions to be equal to infinity
% and return to the calling program
if (sum(params == NaN * ones(k,1)) > 0)
    wtd_moments = Inf;
    return
    disp('Error: did not return');
    break
end


% Error check
if params == NaN * ones(k,1)
   error('params == NaN');
end


% =====================================================================
% Stage 2: Development
% =====================================================================

% Deterministic component of profit from developing
pi0_d = ...
   gamma_tote * tot_e_t_vec + ...
   gamma_totd * tot_d_t_vec + ...
   gamma_bid_bin * bid_bin_t_vec + ...
   gamma_cd_bin * cd_bin_t_vec + ...
   gamma_oilp_bin * oilp_bin_t_vec + ...
   gamma_constant * ones(n_tuples,1);
   

% Initialize Vce
Vce = zeros(n_tuples,Td_m);   

% --------------------------------------
% Form estimates of Vce after lease term
% --------------------------------------

% Vce_I:
% Estimate Vce by inversion (method I)
if (strcmp(Vce_estimator, 'I'))
  Vce_I = ...
    inv(eye(n_tuples) - beta*Me_emp(:,:,Td_m)) * ...
    Me_emp(:,:,Td_m) * gd_emp(:,Td_m) * sigma_eps;
  Vce(:,Td_m) = Vce_I;
end


% Vce_G:
% Estimate Vce by solving for a fixed point after substituting in
% functional form assumptions for expected profits from developing 
% conditional on exploration and for the policy function g_d
% (method G, for the policy fxn)
if (strcmp(Vce_estimator, 'G'))

    % Error check
    disp('Vce_fixed_point');

    % Initial guess for Vce_G comes from 
    % E[profits_if_d] = beta * Vce + sigma_eps
    Vce_G2 = (pi0_d - (sigma_eps * ones(n_tuples,1))) / beta;

    % Initialize value for distance
    distance_VceG = ini_diff_fp;

    % Iterate until convergence
    while distance_VceG > tol_fp 

       % Set function 1 to the previous value of function 2 
       Vce_G1 = Vce_G2;

       % Solve for Vce_G2 as function of Vce_G1
       Vce_G2 = ...
      	  Me_emp(:,:,Td_m) * ...
            (beta * Vce_G1 + ...
             (exp(-(beta * Vce_G1 - pi0_d)/sigma_eps) * sigma_eps));

       % Determine the distance between the two functions 
       % using the sup norm
       distance_VceG = max(abs(Vce_G1 - Vce_G2));
       distance_VceG

       % If distance is NaN, then 
       % set the weighted moment conditions to be equal to infinity
       % and return to the calling program
       if (distance_VceG == NaN)
         wtd_moments = Inf;
         return
         disp('Error: did not return');
         break
       end

    end

    % Determine the continuation value function 
    Vce_G 		= Vce_G2;
    Vce(:,Td_m)		= Vce_G;


end    

% ---------------------------------------
% Form estimates of Vce before lease term
% ---------------------------------------
    
% Loop for each year of lease
for t = Td-1:-1:0

	% Determine matlab index
	t_m = t + 1;

	% Fill in Vce
        Vce(:,t_m) = ...
      	  Me_emp(:,:,t_m) * ...
            (beta * Vce(:,t_m+1) + ...
             (exp(-(beta * Vce(:,t_m+1) - pi0_d)/sigma_eps) * sigma_eps));
	
end	



% --------------------
% Form estimates of gd
% --------------------

% Initialize gd 
gd = zeros(n_tuples,Td_m);


% Loop for each year of lease
for t = 0:Td 

	% Determine matlab index
	t_m = t + 1;

	% Fill in gd for before end of lease term 
        gd(:,t_m) = exp(-(beta * Vce(:,t_m) - pi0_d)/sigma_eps);
	
end	

% Make sure gd is a probability
gd = min(gd,1);



% -------------
% Ex ante value
% -------------
% Ex ante expected value of explored but undeveloped tract,
% where expectations are taken over epsilon.
% Use chosen estimator of Vce (and gd). 
EVe = beta * Vce(:,1) + gd(:,1) * sigma_eps;


% ---------------------------
% Expected conditional profit 
% ---------------------------
% Expected profit from developing conditional on developing, using
% chosen estimator of Vce;
Eprofitifd = beta * Vce + sigma_eps;


 



% =====================================================================
% Stage 1: Exploration 
% =====================================================================

% Estimate deterministic component of profit from exploration
% using chosen estimator of Vce (and gd)
pi0_e = beta * Vce(:,1:Te_m) + ...
        gd(:,1:Te_m) * sigma_eps + ...
        ce_factor * ...
	   (repmat(ce_bin_t_vec,1,Te_m)+ones(n_tuples,Te_m));



% ---------------------
% Form estimates of Vcn
% ---------------------

% Initialize Vcn
Vcn = zeros(n_tuples,Te_m);

% Loop for each year of lease (except the last year,
% when the continuation value is zero)
for t = Te-1:-1:0

	% Determine matlab index
	t_m = t + 1;

	% Fill in Vcn
        Vcn(:,t_m) = ...
      	  Mn_emp(:,:,t_m) * ...
            (beta * Vcn(:,t_m+1) + ...
             (exp(-(beta * Vcn(:,t_m+1) - pi0_e(:,t_m+1))/sigma_mu) ...
	      * sigma_mu));
	
end	


% --------------------
% Form estimates of ge
% --------------------

% Initialize ge 
ge = zeros(n_tuples,Te_m);

% Loop for each year of lease
for t = 0:Te 

	% Determine matlab index
	t_m = t + 1;

	% Fill in ge 
        ge(:,t_m) = exp(-(beta * Vcn(:,t_m) - pi0_e(:,t_m))/sigma_mu);
	
end	

% Make sure ge is a probability
ge = min(ge,1);




% -------------
% Ex ante value
% -------------
% Ex ante expected value of unexplored tract,
% where expectations are taken over mu. 
% Use chosen estimator of Vcn (and ge).
EVn = beta * Vcn(:,1) + ge(:,1) * sigma_mu;




% =====================================================================
% Form moments 
% =====================================================================

% Match predicted to actual explorers (in %)
%    mom_e = ge(tuple_j, params) - I_e_j(tuple_j) * I_not_e_j(tuple_j)
mom_e_vec = zeros(n_tuples,1);
for t=0:Te
  t_m = t + 1;	
  mom_e_vec = mom_e_vec + ...
     tot_not_e_vec(:,t_m) .* (ge(:,t_m) - ge_emp(:,t_m));
end	
mom_e = 1/n * sum(mom_e_vec);

% Match predicted to actual developers (in %)
%    mom_d = gd(tuple_j, params) - I_d_j(tuple_j) * I_e_not_d_j(tuple_j) 
mom_d_vec = zeros(n_tuples,1);
for t=0:Td
  t_m = t + 1;	
  mom_d_vec = mom_d_vec + ...
     tot_e_not_d_vec(:,t_m) .* (gd(:,t_m) - gd_emp(:,t_m));
end
mom_d = 1/n * sum(mom_d_vec);


% Match predicted to actual profits conditional on developing
%    mom_prof = ...
%      E[profits(tuple_j, params) | I_d_j = 1] - profits_j * I_d_j(tuple_j)
mom_prof_vec = zeros(n_tuples,1);
for t=0:Td
   t_m = t + 1;
   mom_prof_vec = mom_prof_vec + ...
     n_d_vec(:,t_m) .* (Eprofitifd(:,t_m) - avg_profitsifd_emp(:,t_m)); 
end
mom_prof = 1/n * sum(mom_prof_vec) / 10;


% Interact mom_e with tot_e_t
mom_e_tote_vec = mom_e_vec .* tot_e_t_vec;
mom_e_tote = 1/n * sum(mom_e_tote_vec);

% Interact mom_e with tot_d_t
mom_e_totd_vec = mom_e_vec .* tot_d_t_vec;
mom_e_totd = 1/n * sum(mom_e_totd_vec);

% Interact mom_e with bid_bin_t
mom_e_bid_bin_vec = mom_e_vec .* bid_bin_t_vec;
mom_e_bid_bin = 1/n * sum(mom_e_bid_bin_vec);

% Interact mom_e with cd_bin_t
mom_e_cdbin_vec = mom_e_vec .* cd_bin_t_vec;
mom_e_cdbin = 1/n * sum(mom_e_cdbin_vec);

% Interact mom_e with oilp_bin_t
mom_e_oilp_bin_vec = mom_e_vec .* oilp_bin_t_vec;
mom_e_oilp_bin = 1/n * sum(mom_e_oilp_bin_vec);


% Interact mom_d with tot_e_t
mom_d_tote_vec = mom_d_vec .* tot_e_t_vec;
mom_d_tote = 1/n * sum(mom_d_tote_vec);

% Interact mom_d with tot_d_t
mom_d_totd_vec = mom_d_vec .* tot_d_t_vec;
mom_d_totd = 1/n * sum(mom_d_totd_vec);



% Form the sample moment vector
moment_vec = [	mom_e;
 		mom_d;	
		mom_prof;
                mom_e_tote;
                mom_e_totd;
                mom_e_bid_bin; 
                mom_e_cdbin;  
                mom_e_oilp_bin; 
                mom_d_tote];
                %mom_d_totd];
moment_vec

% Calculate the weighted moment conditions (1 X 1)
wtd_moments = moment_vec' * inv(weight_matrix) * moment_vec;
wtd_moments

% If any of the moments are NaN, then 
% set the weighted moment conditions to be equal to infinity
% and return to the calling program
if (sum(moment_vec == NaN * ones(length(moment_vec),1)) > 0)
    wtd_moments = Inf;
    return
    disp('Error: did not return');
    break
end







